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1.  Introduction 

The  minimum  requirement  for  the  numerical  evaluation  of  a  flowfield  is 
an  adequate  combination  of  a  physical  description,  discrete  grid  distribu¬ 
tion  and  procedural  algorithm.  All  three  must  allow  for  events  throughout 
the  domain  of  interest.  For  realistic,  complex  flows  this  imposes  varying 
demands  on  each  since  the  disturbance  field  is  generally  nonuniform  and 
frequently  exhibits  local  concentrations  where  there  are  major  changes  in 
state  properties.  It  is  natural  to  attempt  adjustment  of  either  the 
physical  description,  the  grid  pattern  and/or  the  algorithm  itself  so  that 
in  some  sense  there  is  a  match  of  the  available  numerical  power  to  those 
smaller  subdomains  and  the  specific  kinds  of  behavior  found  there.  Such 
considerations  also  appear  on  a  global  scale  when  choices  are  made  of  in- 
viscid  or  viscous  descriptors,  or  of  grid  generation  techniques,  and  of 
methods  that  are  specifically  qualified  for  hyperbolic  or  other  systems. 

The  CFD  challenge  is  to  devise  ways  in  which  appropriate  subdivisions 
of  a  global  domain  may  be  selected  and  if  need  be  altered,  and  how  the 
resulting  regions  may  be  linked  across  their  common  interfaces  during  a 
computation.  Efficient  and  robust  algorithms  must  be  expected  to  reduce 
both  the  computational  speed  and  the  memory  requirements  in  proportion  to 
the  reduction  in  grid  nodes  and  the  simplicity  of  the  description.  Ideally, 
coarse  grids  and  a  minimum  of  physics  are  favored.  Our  focus  is  on  pro¬ 
cedures  which  will  isolate  and  couple  different  portions  of  a  complex  domain 
such  that  in  each  subdomain  the  evaluation  will  disclose  the  details  of  the 
scale  events  to  a  sufficient  and  relevant  accuracy. 

2.  Research  Objectives  and  Tasks 

The  specific  purpose  of  the  research  is  to  develop  procedures  which 
will  permit  subdivision  of  the  field  and  description,  focus  the  computation 
on  those  portions  which  require  a  more  precise  evaluation,  and  control  the 
coupling  between  the  subdomains  so  as  to  maintain  accuracy.  Some  problems 
have  clear  regions  in  which  the  dominant  physics  allows  pre-selection  of 
subdomains.  We  refer  to  this  as  a  non-adaptive  procedure  in  contrast  to 
an  adaptive  approach  in  which  the  evolving  events  guide  and  change  the 


The  program's  initial  effort  began  with  a  determination  of  interface 
constraints  and  necessary  adjustments  of  the  background  Ni  (finite  volume, 
multiple  grid)  method.  Thereafter,  threshold  and  decision  parameters  were 
developed  for  adaptation,  and  alternate  algorithms  have  been  considered  in 
order  to  generalize  the  allowable  mesh  topology,  the  specific  ba^ic 
integrator,  and  the  concept  of  embedding  on  zonal  subdivision  of  an  overall 
domain.  Example  evaluations  have  emphasized  cascade  and  airfoil  flowfields. 

A  brief  summary  of  the  last  year's  efforts  is  described  in  Sections  3,  4 
and  5,  and  some  details  for  several  tasks  are  included  in  attached  Appendices 

3.  Nonadaptive  Embedded  Subdomains 


In  last  year's  annual  report  research  was  reported  on  an  algorithm 
which  it  is  hoped  will  combine  the  best  features  of  node-based  methods 
such  as  Ni's  and  multistage  methods  such  as  Jameson's.  Work  has  continued 
throughout  this  year,  and  good  progress  has  been  made  in  developing  suitable 
damping  and  distribution  operators.  A  detailed  energy  analysis  was  completed 
which  identified  requirements  for  these  operators,  and  eventually  led  to 
suitable  formulations  for  them.  We  refer  the  reader  to  Appendix  A  of  last 
year's  report  for  the  basic  outline  of  the  algorithm.  We  have  not  included 
details  of  the  further  developments  of  this  year,  but  instead  show  some 
results.  Figure  3-1  shows  a  grid  and  Mach  number  plot  for  a  transonic  flow 
over  a  circular  arc  bump  in  a  channel.  The  grid  has  an  embedded  region 
around  the  body.  The  computed  results  are  in  reasonable  agreement  with 
those  of  Dannenhoffer  and  Baron  using  Ni's  method.  Figure  3-2  shows 
supersonic  flow  past  a  4%  bump  in  a  channel.  The  Mach  number  plots 
compare  favorably  with  those  of  Ni  for  this  problem.  The  basic  algorithm 
seems  to  now  be  working.  As  yet  we  have  not  attempted  a  multigrid 
version  which  would  significantly  improve  the  rate  of  convergence. 

However,  we  feel  that  it  is  now  more  important  to  attempt  some  calculations 
with  this  algorithm  for  which  previous  algorithms  had  difficulty.  As 
discussed  last  year,  this  algorithm  has  been  formulated  without  any 
assumption  of  regularity  of  the  mesh  or  the  number  of  dimensions.  We 
have  just  started  to  try  a  calculation  for  a  cylinder  using  the  mesh  shown 
in  Figure  3-3.  It  should  be  noted  that  this  is  a  very  severe  test  of  any 
algorithm  for  the  Euler  equations,  and  we  might  encounter  significant 
problems.  However,  we  think  it  is  important  to  see  how  irregular  a  mesh 
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can  really  be  handled  in  order  to  assess  the  feasibility  of  a  general  grid 
approach.  The  current  algorithm  is  written  to  be  suitable  for  vector  or 
parallel  processors. 

Approach  which  has  less  technical  risk  and  may  prove  to  be  more  prac¬ 
tical  than  the  above  is  the  zonal  approach.  The  flowfield  is  divided  into 
different  zones,  and  within  each  zone  a  standard  algorithm  of  some  type  may 
be  employed.  A  zone  could  have  an  Euler  method  on  a  regular  grid  and  an 
adaptive  method  (within  that  zone) ,  a  Navier-Stokes  method,  or  many  other 
of  a  variety  of  suitable  approaches.  The  important  feature  is  that  the 
zones  communicate  through  the  boundary  data.  Some  thought  has  been  given 
this  year  as  to  how  such  a  zonal  method  should  be  organized  to  be  suitable 
for  multiprocessor  architectures.  A  general  strategy  has  been  conceived, 
but  not  yet  developed  in  detail  nor  tested.  This  will  be  done  during  the 
first  six  months  of  the  next  grant  year. 

As  part  of  the  zonal  approach,  an  idea  for  a  more  efficient  Navier-Stokes 
algorithm  has  been  conceived.  It  is  based  upon  the  observation  that  Navner- 
Stokes  approaches  generally  require  very  fine  grid  spacing  normal  to  the 
body  when  compared  to  the  other  two  directions.  This  leads  to  a  stiff  set 
of  discrete  equations,  which  then  requires  an  implicit  method  for  solution. 
The  standard  approach  is  to  use  a  fully  implicit  method  which  requires 
approximate  factorization  (e.g.  Beam-Warming) .  A  more  attractive  approach 
would  be  to  use  an  implicit  method  only  for  the  differences  in  the  body 
normal  direction,  and  explicit  approaches  for  the  other  two  directions. 

This  algorithm  would  require  a  block  tridiagonal  inversion,  but  no 
factorization.  A  stability  analysis  is  given  in  Appendix  A  for  a  model 
problem,  and  the  results  look  very  promising.  We  will  pursue  this  further 
in  the  next  grant  year. 

4.  Adaptive  Embedded  Subdomains 

The  goal  of  the  adaptive  concept  segment  of  the  effort  has  been  to 
minimize  computational  work  by  appropriately  locating,  controlling  and 
managing  the  embedded  fine  grid  regions.  The  details  involve  procedures 
which  recognize  flow  features  as  they  develop,  decide  on  adequate  thresholds 
for  adaptation  to  be  carried  out,  and  modify  the  discreteness  formulations 
in  order  to  proceed  with  irregular  and  multiple  embedded  grids. 


During  the  first  part  of  the  last  year  the  initial  applications  to  two- 
dimensional  airfoils  were  completed.  Those  calculations  included  a  sensi¬ 
tivity  study  of  the  influence  of  different  refinement  parameters  and  their 
rates  of  change  as  the  basis  for  adaptation  decisions.  Primitive  parameters 
such  as  density,  pressure,  velocity  and  entropy  were  used  as  threshold 
selectors  on  the  basis  of  their  first  and  second  difference  distributions 
throughout  the  domain.  Field  evaluations  for  NACA  0012  and  RAE  2822  lifting 
airfoils  at  transonic  speeds  indicated  density  first  differences  to  be  most 
efficient  and  therefore  preferable  in  defining  subdomains  with  dominant 
activity.  This  extension  of  the  work  to  two-dimensional  fields  was  reported 
on  in  AIAA  Paper  85-0484  and  is  included  here  as  Appendix  C. 

The  effort  has  continued  with  a  further  development  of  the  algorithm 
so  as  to  introduce  a  capability  for  additional  self  evaluation.  Previously 
consideration  was  given  to  the  type  of  feature  expected  and  an  arbitrary 
but  fixed  number  of  adaptive  levels.  A  new  strategy  scheme  has  been 
developed  which  allows  for  automatic  selection  of  the  number  of  embedded 
regions,  in  the  multiple  sense  of  sub-embedding  as  well  as  for  distinct 
features .  A  rule  based  system  monitors  the  progress  and  measures  the  CFD 
needs  against  such  criteria  as  the  changes  in  a  global  parameter  of  interest 
A  number  of  test  cases  have  been  completed  for  airfoil  configurations  that 
have  been  considered  by  the  AGARD  Working  Group  07.  Others  are  planned. 

Some  preliminary  results  are  shown  in  Figures  4-1  through  4-4.  The 
initial  and  final  grids  (convergence  decided  five  adaptations)  for  the 
NACA  0012  airfoil  at  M  =  0.85  and  1.0°  incidence  [AGARD  WG  07  Test  Case  2] 
appear  in  Figure  4-1.  Those  finest  embedded  regions  which  extend  outward 
from  the  unper  and  lower  surfaces  indicate  shock  locations.  The  adaptation 
process  includes  an  allowance  for  growth  of  all  embedded  regions  during  the 
final  stages  of  an  evaluation.  Note  that  a  non-adapted  grid  which  is 
globally  refined  in  order  to  attain  the  same  feature  resolution  would 
contain  a  factor  of  25  times  more  nodes. 

The  surface  Mach  number  distributions  corresponding  to  the  grids  in 
Figure  4-1  are  shown  in  Figure  4-2.  The  Mach  number  contours  in  the  field 
and  convergence  histories  for  lift  and  drag  coefficients  are  ^n  Figure  4-3. 
Lastly,  results  for  M=0.95  are  shown  in  Figure  4-4.  Supersonic  (M=1.2) 


b)  Fifth  (final)  grid  adaptation 


Fig.  4-1  Solution  grids 
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cases  with  detached  bow  and  trailing  edge  fishtail  shocks  also  have  been 
completed.  The  algorithm  is  robust  and  accurate  with  virtually  no  dis¬ 
cernible  effect  due  to  the  sub-grid  interfaces  (see  Figure  4-2).  Generally, 
an  order  of  magnitude  reduction  in  computing  time  is  achieved  by  introducing 
adaptation  in  place  of  global  refinement. 

A  vectorized  version  of  the  adaptive  code  is  currently  being  developed 
and  tested.  The  Cyber  205  at  the  NASA  Langley  Research  Center  is  being 
used  for  that  purpose.  The  intent  is  to  complete  both  globally  refined 
and  adaptive  solutions  so  as  to  provide  precise  comparisons  as  well  as  a 
demonstration  of  the  suitability  for  vectorization.  In  that  regard  the 
data  structure  associated  with  the  adaptation  procedures  appears  to  reduce 
the  benefits  of  vectorization  only  slightly.  A  paper  is  in  preparation  to 
describe  some  of  the  above  findings. 

All  of  the  above  has  continued  the  use  of  the  Ni  method  as  the  basic 
integrator  for  the  algorithm.  A  separate  study  was  undertaken  to  consider 
the  embedding  technique  with  an  alternative  finite  volume  scheme  [Jameson, 
Schmidt  and  Turkel,  AIAA  81-1259] .  The  motivation  is  a  more  efficient  time 
marching  procedure,  but  the  focus  of  the  study  is  on  the  interface  constraints 
and  the  achievement  of  conservative,  accurate  and  smooth  solutions  in  the 
presence  of  abrupt  mesh  scale  changes.  Figure  4-5  shows  grids  and  Mach 
number  contours  for  coarse,  embedded  and  globally  fine  cases  of  transonic 
flow  (M = 0.675)  past  a  typical  channel  bump  (10%)  configuration.  To  date 
analysis  reveals  that  second  order  accurate,  conservative  flux  formulations 
are  not  possible  at  interfaces .  Smoothing  also  implies  nonconservative 
contributions.  Three  classes  of  interface  formulations  are  suggested: 
zeroth  order  accurate  and  fully  conservative,  first  order  accurate  for 
conservative  fluxes  and  nonconservative  smoothing,  and  second  order 
accurate  but  fully  nonconservative.  The  Mach  number  contours  in  Figure 
4-5  illustrate  the  advantages  of  global  refinement.  The  interface  effects 
are  also  apparent  along  the  superimposed  interface  locus  for  the  embedded 
case.  This  work  with  the  Jameson  scheme  is  continuing. 

5.  Jerusalem  Meeting 

A  workshop  was  held  in  Jerusalem,  Israel  in  December  1984  entitled  "The 
Impact  of  Supercomputers  on  the  Next  Decade  of  CFD."  The  participation  of 
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b)  Global  fine  grid  (64  x 16)  and  corresponding  Mach  number  contours  (AM 
Min  =  0.37,  Max  =  1.64 

Fig.  4-5  Transonic  channel  flow,  M  =  0.675,  10%  thick  bump 


Doubly  embedded  grid,  base  grid  32  x  8 


Fig.  4-5  (concluded) 
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r  Professor  Murman,  one  of  the  organizers  of  the  workshop,  was  sponsored  by 

b  this  grant.  Professors  Murman  and  Abarbanel  have  prepared  a  summary  report 

of  the  workshop  which  is  attached  as  Appendix  C.  A  number  of  important 
issues  relative  to  this  grant  were  brought  out  at  the  workshop  and  are 
recorded  in  Appendix  C.  Perhaps  the  most  important  are  the  observations 
*  that  current  algorithms  cannot  simply  be  scaled  up  to  256  megaword  machines. 

They  then  will  converge  so  slowly  that  the  procedure  will  be  impractical. 
Embedded  and  adaptive  mesh  approaches  look  very  favorable  as  a  way  to 
;  proceed. 
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where  X  =  At/Ax  and  X  =  At/Ay  . 

x  y 

Define  z  =  X  sin  0  +  A  sin  $  ;  the  amplification  factor,  G  =  Gn+^/un  then 
x  y 

becomes : 


G  =  1  +  iz  -  z2  -  iz3  =  (1  -  z2) (i  +  iz) 


(1.9) 


and  therefore 


I G 1 2  =  (1  -  zVd  +  Z2)  =  1  -  Z2  -  Z4  +  Z6 


(1.10) 

The  Von  Neumann  stability  requirement  | G |  -  1,  then  leads  to  the  following 
inequality: 


z4  -  z2  -  1  ^  0  , 


(1.11) 


leading  to: 

)1/2 

~  1.27202  (1.12) 

j 

Recall  now  that  z  =  X  sin  0  +  X  sin  d>,  and  requirement  (1.12)  takes 

x  y 

the  form 


^.+At=X  +X  5  1.2720 
Ax  Ay  x  y 

or,  upon  rearranging 
At  <  1.2720 

Ax  +  Ay 


(1.13) 


(1.14) 


The  stability  condition,  (1.14)  ,  clearly  shows  that  if  in  the  calculation 
we  choose  the  time  step  accordingly,  i.e. 


.  <  min 

At  -  .  , 
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1.2720 


•  ay-i,k 


S.k +  Ayj.k 


(1.15) 


we  may  get  very  severe  restriction  if  the  mesh  is  highly  stretched  so  that 

somewhere  in  the  field  (Ay/ Ax) .  ,  say,  is  very  small. 

3 

To  overcome  this  difficulty  without  paying  an  inordinately  high  price 
of  two  matrix  inversions  per  stage  we  propose  a  semi-implicit  algorithm. 


2.  Semi-Implicit  Scheme 

Using  the  previous  notation,  the  semi-implicit  3-stage  Runga-Kutta 
algorithm  takes  the  form: 
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Again  "telescoping,"  and  joined  directly  to  Fourier  space,  the  amplification 
factor  is  found  to  be: 


G  = 
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where  £  =  sin  0,  r)  =  sin  ;  i.e.  -1  -  £  r)  -  1.  The  Von  Neumann  stability 

requirement  leads  to  the  following  inequality: 
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where  x  =  X  5  and  y  =  Xr|;  i .  e .  -  00  <  <  00  . 

x  y 

We  carry  out  the  investigation  of  (2.7)  by  first  considering  the  case 
of  y2  >  x2.  Divide  (2.7)  by  y2  -  x2  >  0  to  get 


(y  +  x)2 
2  2 


2  2 
x  y 


y  -  x 

The  worst  case  is  that  of  xy  <  0.  Let  x->--x,  y  + 
xy  >  0  and  the  inequality: 


-  y4  ^  0  (2.8) 

y,  then  the  new  x,y  satisfy 
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y  -  x 

2  2 

Note  that  in  (2.9)  it  does  not  matter  (with  y  >  x  and  xy  >  0)  whether  x 
and  y  are  both  negative  or  positive.  We  shall,  therefore,  consider  the 
x  >  0,  y  >  0  case.  We  rewrite  (2.9)  slightly: 


(2.9) 
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y  -  x  ,  2  4 

- -  +  x  -  x 
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5  -  1  +  X2  -  X4 


(2.10) 


The  inequality  (2.10)  is  always  satisfied.  Therefore,  we  have  shown  that 

2  2 

(2.7)  is  satisfied  when  y  >  x  for  all  -  00  <  x,y  <  00  .  Notice  also, 

directly  from  (2.7)  that  when  y  =  -x  the  inequality  is  satisfied  (if  y  = -x 

2  2 

we  have  the  equality).  Therefore,  (2.7)  is  satisfied  for  all  y >  x  , 


Next  we  consider  the  case  x  >  y  .  We  see  from  (2.7)  that  now  the  worst 

case  is  that  of  xy  >  0.  Again  x  and  y  can  both  be  either  negative  or  positive 

without  changing  the  left  hand  side  of  (2.7) .  We  shall  use  x  >  y;  x  >  0, 

2  2 

y  >  0.  Divide  (2.7)  x  -  y  >0: 


(x+  y) 
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.  2  2.  2  4  22  4  <- 

(x  -  y  )  +  y  +  2xy  +  x  +  xy  +  y  5  0 


(2.11) 


Since  x  >  y,  take  y  =  8x,  (0  <  8  <  1) .  Inequality  (2.11)  becomes: 
~  ¥4  +  <262  +  26  -  Dx2  +  (64  +  82  +  l)x4  5  0 


or,  with  w  =  x  , 


/(WjB)  =  (84  +  62  +  1) w^  +  (28^  +  26-Dw  -  <0 


(2.12) 


The  smallest  w,  for  all  0  <  8  <  1/  that  solves  (w.B)  =  0  is  w  =  4/3  for  8 =  1/2 
This  is  the  stability  criterion.  Thus  we  have 

Aj£|  5  (4/3) 1/2 

and  so  the  stability  condition  becomes  At/Ax  =  X  <  2//3,  or 


At  <  1.1547 


mi"  Ax.  J 

D»k  3,kj 


(2.13) 


Notice  that  now  we  have  a  conditional  stability  that  does  not  depend  on 
the  mesh  stretching.  This  condition  is  still  superior  to  (1.15)  even  when 
Ax^  k  =  Ay.,  .  everywhere,  because  then  (1.15)  yields  At  -  0.636  Ax_.  . 
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ABSTRACT 

An  adapted  grid  technique  for  accurate 
and  efficient  solutions  of  the  steady,  two- 
dimensional  Euler  equations  is  presented. 
Embedded  meshes,  which  are  coupled  to  a 
fixed  global  mesh  via  a  multiple-grid 
technique,  are  employed.  The  effect  of 
choosing  various  adaptation  strategies  is 
examined.  Solutions  are  presented  for 
isolated  lifting  airfoils  in  subcritical  as 
well  as  transonic  flow.  The  supersonic  flow 
over  a  circular-arc  cascade  with  a  complex 
shock  structure  is  also  presented.  In  all 
cases,  the  adapted  solution  achieved  the 
same  accuracy  as  global  refinement,  but 
required  a  factor  of  between  5  and  7  less 
computer  time  and  between  3  and  8  less 
storage. 

INTRODUCTION 

The  numerical  solution  of  the  two- 
dimensional  Euler  equations  is  typically 
carried  out  by  discretizing  the  flow  field 
and  then  solving  a  coupled  set  of  equations 
for  each  node.  This  technique  suffers  from 
the  introduction  of  errors  due  to  the 
discrete  approximation  of  the  continuous 
governing  equations.  These  errors  are 
related  to  the  local  mesh  spacing  and  the 
local  solution  smoothness. 

To  obtain  more  accurate  results,  one 
can  decrease  the  mesh  spacing.  But  since 
the  domain  size  is  fixed,  this  results  in 
more  required  grid  points,  with  the 
associated  Increase  in  computational  work 
and  required  storage.  It  has  been  noted  by 
many  prior  research ers  C13  that  for  the 
Euler  equations  applied  to  problems  of 
Interest,  much  of  the  flow  field  is  smooth, 
resulting  in  small  errors.  There  are 
however  Isolated  areas  in  the  flow  field 
where  the  solution  is  not  smooth,  resulting 
in  appreciably  larger  errors.  Thus,  it  is 
an  efficient  strategy  to  refine  the  mesh 
only  where  necessary. 
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Unfortunately,  one  does  not  in  general 
know  £  priori  where  these  unsmooth  regions 
are,  and  thus  it  often  takes  a  few  trials 
and  errors  to  put  the  dense  mesh  where  it  is 
really  required.  To  circumvent  this 
problem,  adaptive  grid  techniques  have  been 
developed  to  search  for  regions  of  high 
errors,  and  then  refine  the  mesh  there. 
Thus,  the  solution  and  the  grid  evolve 
simultaneously. 

There  are  currently  two  methods  for 
adapting  grids.  The  first  (and  most 
popular)  is  the  grid  point  redistribution 
scheme  which  has  been  reviewed  extensively 
by  Thompson  C13.  In  this  technique,  there 
is  a  fixed  number  of  grid  points  which  move 
around,  congregating  in  regions  of  interest. 
It  is  reasonably  easy  to  program,  since  at 
every  iteration  (or  at  least  every  few 
iterations),  some  statistic  of  the  flow 
field  is  sensed  (for  example  entropy 
gradient),  and  the  grid  points  move  either 
by  an  attraction  model  or  by  using  a  Poisson 
equation  grid  generator.  Redistribution 
however  has  a  few  drawbacks.  First,  since 
there  la  a  fixed  number  of  grid  points, 
moving  them  toward  flow  features  deprives 
other  flow  -regions  of  adequate  resolution. 
Thus,  one  could  say  that  the  solution 
becomes  both  good  uniformly  and  at  the  same 
time  uniformly  bad.  Second,  due  to 
topological  restrictions,  grid  lines  which 
properly  congregate  near  a  shock  at  an 
airfoil  surface  must  propagate  out  into  the 
flow  field,  often  to  the  far  field  boundary 
where  the  increased  resolution  is 
unnecessary.  Thompson  also  points  out  that 
excessive  grid  skewing  and  stretching  can 
occur,  producing  Inherent  inaccuracies. 

The  second  grid  adaptation  technique 
involves  locally  embedding  sub-grids  in 
regions  of  interest.  This  results  in  an 
Increase  in  the  number  of  grid  points,  and 
hence  an  increase  in  required  computational 
resources.  Also,  artificial  internal 
boundaries  are  created  in  the  flow  field; 
care  must  be  taken  to  ensure  that  the 
solution  will  smoothly  transfer  across  such 
an  interface,  while  maintaining  conservation 
and  stability.  Since  global  grid  points  are 
never  moved,  the  accuracy  is  not  reduced 
away  from  flow  features  as  in  redistribution 
methods. 


Berger  and  Jameson  C23  have  used  local 
grid  embedding  for  the  two-dimensional  Euler 
equations.  In  their  approach,  the  flow 
features  are  detected  and  "best-fit" 
rectangles  are  superimposed  over  the  global 
grid.  This  results  in  a  mesh  interface 
across  which  it  is  very  difficult  to  enforce 
conservation.  In  addition,  determination  of 
the  best-fit  rectangle  requires  complex 
clustering  algorithms. 

An  alternative  embedded  mesh  procedure 
has  been  proposed  by  Dannenhoffer  and  Baron 
C31.  In  this  approach,  irregularly  shaped 
embedded  regions  which  are  topologically 
connected  to  the  global  grid  are  used.  This 
results  in  both  a  simpler  embedded  mesh 
interface  as  well  as  a  simpler  embedding 
scheme  which  doesn't  require  clustering. 


The  use  of  Irregularly  shaped  embedded 
regions  introduces  a  bookkeeping  problem 
which  must  be  addressed.  Since  Ni's  scheme 
is  cell-based,  each  cell  can  be  integrated 
independently  which  is  an  important 
consideration  in  the  choice  of  this  method. 
Cells  communicate  with  each  other  only 
through  the  variables  at  the  shared  nodes 
from  the  previous  explicit  pseudo-time  step. 
A  data  structure  which  is  based  on  that 
property  was  first  introduced  by  Usab  and 
Murman  C53,  a  variation  of  which  is  used  in 
the  current  program. 

GOVERNING  EQUATIONS 

The  unsteady,  two-dimensional  Euler 
equations  in  conservation  law  form  are  given 
by 


In  reference  C33,  the  authors 
discussed  their  technique  for  the 
one-dimensional  Euler  equation  and 
two-dimensional  Burgers  equation.  The 
extension  of  their  technique  to  the 
two-dimensional  Euler  equations  is  presented 
here. 


This  paper  first  discusses  the  overall 
adaptation  strategy.  The  governing 
equations  are  presented  and  the  embedded 
mesh  algorithm  is  described.  Included  in 
that  section  are  discussions  of  the  embedded 
mesh  interface  as  well  as  the  adapted 
smoothing  algorithm  employed.  Appropriate 
choices  for  a  refinement  parameter  are 
discussed  in  the  next  section.  Computed 
results  for  the  Euler  flow  over  airfoils  and 
a  model  problem  with  complicated  shock 
topology  complete  the  discussion. 
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In  equation  (1),  p  is  the  density,  u  and  v 
are  the  velocity  components  in  the  x-  and  y- 
directions,  E  is  the  total  internal  energy 
per  unit  volume,  p  is  the  pressure,  h  is 
the  total  enthalpy,  and  y  is  the  ratio  of 
specific  heats. 


ADAPTIVE  SOLUTION  ALGORITHM 


riME  MARCHING  PROCEDURE 


The  basic  approach  of  the  adapted 
solution  procedure  described  here  is  to  use 
a  fixed  global  grid  and  then  embed 
irregularly  shaped  grids  where  necessary. 
This  is  accomplished  by  solving  initially 
only  on  the  global  grid.  When 
quasi -convergence  is  reached,  a  refinement 
parameter  is  computed  at  each  node  and  any 
cells  which  are  connected  to  nodes  where  the 
refinement  parameter  is  above  some  threshold 
are  divided.  The  program  then  reintegrates 
on  the  new  grid  (global  and  embedded, 
coupled  by  a  multiple-grid  scheme).  After  a 
prespecified  number  of  adaptations 
(typically  2),  the  program  iterates  to  the 
final  converged  solution. 

The  integration  of  the  governing 
equation  is  performed  using  Ni's  multiple- 
grid  algorithm  C43.  This  scheme  is  composed 
of  two  parts  —  a  finite-volume  form  of  a 
standard  single-step,  Lax-Wendroff 
integration  applied  on  a  fine  mesh,  and  a 
coarse  grid  accelerator  which  operates  on 
residuals  transported  from  the  fine  mesh 
solver.  In  both  parts,  a  “change”  is 
computed  at  the  center  of  each  cell  and  then 
transferred  to  the  adjacent  nodes  by  means 
of  ‘distribution  formulae”. 


Basic  Scheme 

Consider  the  fine  mesh  cell  shown  in 
figure  1.  To  calculate  the  "change"  in  the 
dependent  variables  at  the  center  of  this 
cell,  the  divergence  theorem  is  applied  to 
the  governing  equations,  giving 

"  *5  (  fw  -  +  )  (2> 

where  F  denotes  the  contravariant  flux 
through  each  given  face,  dt  denotes  the 
pseudo-time  step,  and  dV  the  cell  volume. 
The  contravariant  fluxes  are  computed  by 
trapezoidal  integration  along  each  cell 
face,  based  on  the  dependent  variables  at 
the  corner  nodes.  Equation  (2)  is  called 
the  cell  flux  balance. 


Figure  1.  Typical  computational  cell. 


The  distribution  formulae  serve  to 
transfer  the  "change"  <AU>  from  the  center 
of  the  cell  to  the  four  corner  nodes.  The 
foraulae  are  derived  froa  a  second-order- 
accurate  Taylor  series  expansion  of  U  (with 
respect  to  tlae) ,  and  are  given  by 

“  i  (  *u  |  s  ^  |  £  AG  )  <3) 

where 


the  characteristic  variables  are  predicted 
correctly  by  the  interior  scheae. 

Eabedded  Mesh  Interface 

The  embedded  grid  is  composed  of  cells 
which  were  formed  by  dividing  a  global  (or 
previously  divided)  cell.  Locally,  this 
appears  exactly  like  the  fine  and  coarse 
grid  cells  used  in  the  multiple-grid 
algorithm. 


AF  -  fff  dU  and  AG  -  AU  (4) 


are  the  unsteady  fluxes  based  upon  the 
Jacobians  of  F  and  C  evaluated  at  the  center 
of  the  cell.  The  first  term  in  equation  (3) 
(AU)  is  the  first-order-change-in-time  for 
the  Taylor  series  expansion  while  the  last 
two  terms  represent  a  second -order-change- 
in-time  which  is  necessary  for  stability. 
The  latter  terms  bias  the  distribution  of 
the  "change"  in  the  windward  direction, 
which  is  similar  to  the  stabilizing  effects 
of  upwind  differencing  £63. 


To  accelerate  solution  convergence,  Ni 
Introduced  a  multiple-grid  algorithm  £43 
which  propagates  the  fine  grid  changes  at 
coarse  grid  speeds.  The  coarse  grid  is 
generated  by  eliminating  every  other  fine 
grid  line.  The  coarse  grid  acceleration  is 
accomplished  by  transporting  the  changes 
previously  calculated  from  the  fine  meshes 
to  the  center  of  a  coarse  grid  cell  and  then 
distributing  that  change  to  the  coarse  grid 
nodes  by  use  of  the  distribution  formulae 
( 3) . 

Boundary  Conditions 

For  the  problems  considered  in  this 
paper,  there  are  two  types  of  boundary 
conditions.  The  first  is  a  solid  wall 
boundary  condition  where  "no  flow  through 
the  surface"  is  enforced.  If  one  views  the 
solid  wall  as  a  streamline,  then  the  effect 
of  a  pseudo-cell  just  Inside  the  body  can  be 
considered.  By  combining  the  effects  of  the 
distribution  formulae  (3)  on  each  side  of 
the  wall,  one  finds  that  the  boundary 
condition  reduces  to  doubling  the  boundary 
changes  predicted  by  the  interior  cell,  and 
cancelling  the  normal  velocity  component. 

The  other  boundary  condition  type  is  a 
far  field  condition.  For  an  airfoil 
calculation,  the  far  field  is  composed  of  a 
uniform  flow  plus  the  effect  of  a  vortex, 
whose  strength  is  set  based  upon  an 
approximate  lift  coefficient.  For  the 
cascade  problem,  the  free  stream  is  assumed 
to  be  simply  a  uniform  flow  both  upstream 
and  downstream. 

The  boundary  conditions  are  applied  at 
far  field  nodes  by  using  a  characteristic 
analysis  in  the  local  streamwise  direction; 
the  characteristic  variables  which  enter  the 
domain  remain  unchanged,  while  for  those 
exiting,  it  is  assumed  that  the  changes  in 


Using  a  multiple-grid  accelerator  to 
couple  the  global  and  embedded  meshes  was 
first  suggested  by  Brandt  £73  and  was 
implemented  by  Brown  for  the  full  potential 
equation  CB3.  With  this  technique,  waves 
propagate  through  the  embedded  regions  at 
coarse  grid  speeds.  Usab  has  shown  that 
this  coupling  results  in  convergence  rates 
which  are  as  fast  as  coarse-grld-alone 
solutions  £93.  This  is  significant  when  one 
considers  the  coi..-  *quences  of  simply 
coupling  global  and  embedded  regions  at  the 
interface  £23.  In  the  latter  technique, 
wave  propagation  is  restricted  to  the 
embedded  (fine)  grid  speed,  resulting  in 
slower  convergence  rates. 

Points  at  the  edge  of  the  embedded 
domain  must  be  carefully  treated  in  order  to 
maintain  global  conservation  and 
computational  stability.  Consider  the 
embedded  mesh  interface  shown  in  figure  2. 
In  the  present  scheme,  nodes  2,  5,  and  9  are 
considered  part  of  the  fine  domain.  Thus 
the  flux  balance  and  distribution  formulae 
can  be  applied  as  usual  to  cells  B,  C,  D, 
and  E.  Although  the  changes  in  the 
dependent  variables  are  computed  at  nodes  2, 
5,  and  9  due  to  cells  B  and  E,  the  changes 
must  not  be  applied  when  operating  on  the 
embedded  mesh  level.  Instead,  they  must  be 
stored  and  applied  only  after  the  (explicit) 
flux  balance  has  been  performed  on  cell  A. 

Since  cell  A  is  a  fine  cell  on  the 
global  grid  level,  the  appropriate 
integration  consists  of  a  flux  balance  and 
distribution.  Each  of  these  steps  must  be 
modified  due  to  the  presence  of  node  5.  The 
flux  balance  now  consists  of  the  sum  of  the 
fluxes  through  five  faces  (1-2,  2-5,  5-9, 
9-8,  and  8-1).  This  can  be  easily 
incorporated  into  the  trapezoidal 
integration.  Since  the  fluxes  through  the 
interface  (for  example  2-5)  cancel  in  the 
adjoining  cells  (for  example  A  and  B),  the 
scheme  maintains  globs')  conservation. 


Figure  2.  Detail  of  embedded  mesh. 
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The  distribution  formulae  given  by 
equation  <3)  remain  valid  for  cell  A. 
However,  the  change  at  the  center  of  the 
cell  must  also  be  distributed  to  node  5. 
This  is  accomplished  by  averaging  the 
distribution  to  nodes  2  and  9,  or  in  general 
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Coarse  cell  BCDE  Is  treated  the  same 
as  any  other  coarse  cell  interior  to  the 
embedded  region.  This  yields  an  apparent 
Inconsistency  at  nodes  2,  5,  and  9  due  to 
the  absence  of  a  coarse  cell  underlying  cell 
A.  At  convergence  however,  the  residuals 
transferred  to  BCDE  do  vanish  as  do  the 
Inconsistencies. 

Adapted  Smoothing 

As  with  most  solutions  of  the  Euler 
equations,  artificial  viscosity  is  needed  to 
damp  out  spurious  oscillatory  solutions  and 
to  capture  shocks.  Different  levels  of 
artificial  viscosity  are  needed  for  these 
two  purposes.  Jameson  CIO]  uses  fourth 
order  smoothing  over  the  whole  flow  field  to 
damp  out  spurious  oscillations,  with  second 
order  smoothing  Blended  in  near  shocks  to 
capture  them.  This  blending  has  the  effect 
of  putting  the  lower  order  smoothing  only 
where  It  is  needed. 


An  alternative  way  of  using  the  proper 
smoothing  only  where  needed  Is  to  use  second 
order  smoothing  globally,  but  to  vary  the 
smoothing  coefficient  so  that  only  the 
required  amount  is  used  locally.  Typically, 
the  coefficient  Is  related  to  the  second 
difference  of  density  or  pressure. 

For  Illustrative  purposes,  the 
following  will  be  developed  in  one 
dimension,  with  Its  extension  to  two 
dimensions  being  straightforward. 

The  smoothing  at  any  point  i  is  given 
by 

smooth^  -  (u^j^Uj+u^^)  <6> 

where  u  is  the  parameter  being  smoothed  and 
p,  is  the  spatially  varying  smoothing 
coefficient  (which  includes  the  appropriate 
geometric  parameters). 


To  vary  the  smoothing  coefficient  from 
node  to  node,  a  second  difference  of  density 
Is  used  in  the  current  program.  This  takes 
the  form 
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Here  n  is  the  level  of  background  smoothing 
needed  to  control  spurious  oscillations  and 
6  Is  the  factor  by  which  the  smoothing  has 
to  be  increased  (above  the  background  level) 
near  shocks.  Typical  values  are  n*0.025  and 
6-50. 


A  desirable  characteristic  of  any 
smoothing  scheme  is  that  it  be  conservative, 
l.e.,  the  sum  of  the  smoothing  terms  over 
the  whole  domain  vanishes.  If  p,  is  a 
constant  in  (6),  the  contributions  of  the 
interior  nodes  in  fact  do  cancel,  leaving 
only  boundary  contributions.  If  however  p± 
varies,  this  is  not  true.  To  circumvent 
this  problem,  one  can  use  an  averaging  of 
the  smoothing  coefficients,  giving 

smoot^  -  +  (p^j+p.)  *  <ui-l”ul) 

(8) 
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Thus,  the  sum  of  the  smoothing  over 
all  i  results  In  the  cancellation  of  the 
contributions  of  all  the  interior  intervals, 
leaving  only  the  boundary  terms,  or 

N-l 

E  (smooth^  *  +  (p2  +p2  )  *  (u2  -u2  ) 

-  *  <UN_i-UN>  <9> 


This  leaves  two  possibilities.  First, 
one  can  accept  the  amount  of  non¬ 
conservation  near  the  boundary  and  allow  the 
creation  of  mass  (or  momentum);  the 
resulting  error  is  not  serious.  In 
practice,  a  more  serious  problem  with  this 
approach  is  that  the  smoothing  operator 
looks  convective  (rather  than  diffusive) 
near  the  boundary.  This  is  not  an 
acceptable  numerical  model  for  a  smoothing 
operator.  The  other  approach  (and  the  one 
adopted  here)  is  simply  not  to  smooth  normal 
to  the  ends  of  the  domain.  As  long  as 
smoothing  is  not  required  there  (and 
numerical  experiments  show  that  it  is  not), 
this  is  a  valid  approach. 

A  spatially  varying  smoothing 
coefficient  can  also  cause  the  smoothing 
term  to  look  convective  within  the  flow 
field.  This  occurs  if  p.  is  allowed  to  vary 
too  rapidly  within  the  field.  It  can  easily 
happen  near  shocks  where  the  second 
difference  of  density  can  change  rapidly 
from  point  to  point.  It  is  reasonable 
therefore  to  require  a  smooth  distribution 
of  p.,  which  requires  smoothing  the 
smoothing  coefficient.  This  can  be 
accomplished  by  calculating  p.  at  every 
point,  and  then  smoothing  the  values  before 
use  in  the  smoothing  step.  Unfortunately, 
this  requires  three  sweeps  through  all  the 
computational  nodes,  which  is  expensive. 
Therefore,  in  the  present  technique,  a 
special  smoothing  coefficient  equation 
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(10) 


is  sdded  to  the  set  of  equations  being 
solved,  where  n  is  the  iteration  count  and 
w.  is  the  desired  smoothing  coefficient 
computed  from  equation  (7).  In  equation 
(10),  a  governs  the  rate  of  convergence  for 
p, ,  which  must  he  faster  than  the  rate  of 
convergence  of  the  Euler  equations.  The 
parameter  0  in  (10)  governs  the  smoothness 
of  This  equation  assumes  that  p4  does 
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not  vary  much  from  the  previous  iteration, 
and  thus  using  the  prior  values  for  the 
diffusive  terms  does  not  cause  significant 
errors.  Typical  values  are  a-0.1  and  0-1. 


COMPUTED  RESULTS 


The  computed  results  are  presented  in 
two  sections.  In  the  first,  various 
refinement  parameters  are  examined,  and 
appropriate  choices  for  the  two-dimensional 
Euler  equations  are  discussed.  The  second 
section  compares  adapted  solutions  with 
solutions  computed  with  global  refinement. 
Here,  global  refinement  refers  to  the 
process  of  subdividing  every  cell  in  the 
domain.  These  comparisons  are  made  with 
respect  to  both  accuracy  and  efficiency. 


Feature  Detection 


Adaptation  requires  that  the  program 
be  able  to  sense  where  and  when  refinement 
is  required.  For  the  Euler  equations  there 
are  many  values  which  might  be  sensed  for 
this  purpose.  For  example,  one  may  use 
first  or  second  differences,  or  even  the 
truncation  error.  Those  operators  can  be 
applied  to  such  primitive  values  as  density 
or  energy,  or  to  derived  quantities  such  as 
the  entropy. 


In  choosing  an  appropriate  refinement 
parameter,  one  must  consider  the  kinds  of 
features  which  are  to  sensed  and  things 
about  them  which  are  different  than  the 
background  flow.  For  example,  if  one  is 
searching  for  shocks,  it  is  clear  that  a 
change  in  density,  pressure,  and  entropy  is 
generated  across  the  shock,  but  that  there 
is  no  change  in  mass.  Similarly,  if  one  was 
looking  for  the  slip  line  behind  an  airfoil, 
it  may  be  useful  to  use  density  or  entropy, 
but  not  static  pressure  since  it  is 
continuous  across  the  slip  line. 


In  addition  to  choosing  which  variable 
to  sense,  one  must  also  select  a  way  of 
neasurlng  how  it  varies.  Typically  first  or 
second  differences  might  be  examined.  Here, 
the  undivided  differences  were  chosen 
Instead  of  the  gradient  or  Laplaclan  of  the 
chosen  variable  since  it  is  important  that 
the  measure  react  to  prior  adaptation.  For 
example,  if  one  used  the  entropy  gradient  as 
a  measure  near  a  shock,  adaptation  would 
seem  to  have  the  wrong  effect,  since  after 
adaptation,  the  gradient  would  Increase, 
indicating  a  need  for  more  adaptation,  etc. 


On  the  other  hand,  if  the  first  difference 
of  density  were  used,  then  the  refinement 
P*raa«ter  would  not  Increase  as  a  result  of 
prior  adaptation. 


Figure  3  shows  the  effect  of  choosing 
different  refinement  parameters  during 
actual  computations.  The  embedded  grids 
were  computed  for  an  RAE  2822  airfoil  at 
Mach  number  0.75  and  angle  of  attack  3.0 
degrees.  This  case  has  a  shock  on  the  upper 
surface  at  about  75  percent  chord.  In  each 
of  the  cases,  there  are  two  sequential 
adaptations,  with  the  threshold  selected 
automatically  as  described  below. 

Figure  3a  shows  the  grid  resulting 
from  using  the  first  difference  of  density 
as  the  refinement  parameter.  In  this  case, 
double  embedding  was  automatically  generated 
around  the  stagnation  point,  the  shock,  and 
the  trailing  edge.  In  addition,  double 
embedding  followed  the  expansion  fan 
generated  at  the  leading  edge.  The 
remainder  of  the  airfoil  surface  is  singly 
embedded,  except  on  the  lower  surface  near 
the  trailing  edge. 


Figure  3b  shows  the  grid  resulting 
from  using  the  first  difference  of  pressure 
as  the  refinement  parameter.  The  grid  is 
very  similar  to  that  in  figure  3a,  except 
that  there  is  less  adaptation  near  the 
trailing  edge  and  in  the  near  wake  region. 
This  is  due  to  the  fact  that  pressure  is 
continuous  across  the  slip  line  emanating 
from  the  trailing  edge,  whereas  the  density 
has  a  jump. 

In  figure  3c,  the  first  difference  of 
velocity  served  as  the  refinement  parameter. 
Here  again,  the  pattern  of  adaptation  is 
very  similar  to  the  previous  two  examples. 
The  main  difference  is  that  the  extent  of 
the  adaptation  behind  the  airfoil  is 
slightly  larger  than  either  of  the  other 
two. 


The  first  difference  of  entropy  is 
used  as  the  refinement  parameter  in  figure 
3d.  It  can  be  seen  that  the  adaptation  very 
faithfully  follows  the  entropy  gradient  set 
up  in  the  wake,  even  out  to  the  far  field 
boundary.  Since  one  must  control  the 
increase  in  required  computational 
resources,  this  results  in  fewer  points 
being  available  for  other  parts  of  the 
field,  with  the  lower  surface  of  the  airfoil 
being  almost  completely  Ignored.  Also,  the 
adaptation  region  around  the  shock  is  rather 
small,  resulting  in  an  appreciable  strength 
shock  passing  through  the  edge  of  the 
embedded  region. 

The  computed  lift  and  drag 
coefficients  for  each  of  these  cases  were 
all  within  1.4  percent  of  each  other, 
indicating  that  the  adaptations  were  equally 
good  from  the  standpoint  of  accuracy.  One 
can  see  from  the  figures  that  the  normalized 
CPU  times  to  convergence  ranged  from  4.5  to 
7.1.  However,  these  times  were  not  directly 


a)  First  difference  of  density 
Nodes  -  1601,  Tine  -  4.5 


e)  Second  difference  of  density 


f)  Second  difference  of  pressure 


c)  First  difference  of  velocity. 
Nodes  «  1641,  Tine  -  4.8 


Figure  3. 

Effect  of  choice  of  refinement  parameter. 
RAE  2822,  Mach-0.75,  a-3.0  deg. 
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related  to  the  number  of  nodes  for  each  slope).  This  corresponds  to  the  decrease  in 
case.  due  to  the  varying  number  of  the  threshold  below  the  value  of  noise  In 
Iterations  required  to  reach  convergence,  figure  4.  If  on  the  other  hand  there  is  a 
From  these  cases  as  well  as  others  which  the  slight  increase  in  the  threshold,  there  are 
authors  have  reviewed,  it  appears  that  the  very  few  nodes  which  would  no  longer  become 
adaptations  which  completely  surround  flow  adapted.  Therefore,  the  'knee*  in  figure  S 


features  (shocks,  leading  and  trailing 
edges)  perform  the  best. 

Figures  3e  through  3h  show  the  effect 
of  using  the  second  difference  of  density, 
pressure,  velocity,  and  entropy 
respectively.  In  each  case,  the  edge  of  the 
embedded  region  proves  to  be  relatively 
ragged.  There  are  voids  in  the  embedded 
region  as  well  as  islands  of  adaptation,  as 
easily  evidenced  in  figure  3f.  This  results 
in  convergence  difficulties  in  the  solver. 
Consequently,  all  of  the  cases  using  second 
differences  resulted  in  unacceptably  long 
computation  times  due  to  such  topological 
problems. 

Based  on  these  findings,  as  well  as  on 
the  accuracies  and  computation  times  for  all 
these  methods,  it  appears  that  the  first 
difference  of  density  is  the  best  refinement 
parameter,  especially  for  transonic  airfoils 
with  shocks,  stagnation  regions,  and  slip 
lines. 

In  addition  to  selecting  an 
appropriate  refinement  parameter,  it  is 
necessary  to  determine  a  threshold; 
adaptation  is  performed  around  any  nodes 
whose  refinement  parameter  exceeds  this 
threshold.  This  is  a  classical 
signal-to-noise  discriminator. 

Figure  4  shows  the  distribution  of 
refinement  parameter  arbitrarily  plotted 
versus  node  number  for  the  airfoil 
calculation  previously  discussed.  The 
refinement  parameter  for  this  case  is  the 
first  difference  of  density,  normalized  by 
the  average  first  difference  over  the  whole 
domain.  There  are  four  levels  of  possible 
threshold  shown. 

Choosing  A  as  the  threshold  results  in 
very  few  cells  being  adapted. 
Alternatively,  choosing  D  results  in  almost 
all  the  cells  being  adapted,  yielding  an 
almost  global  embedding,  which  is  a  very 
ineffective  strategy.  It  appears  from 
figure  4  that  an  appropriate  threshold  value 
would  be  either  B  or  C. 

This  data  is  presented  again  in  figure 
5  which  is  the  cumulative  distribution 
function  corresponding  to  figure  4.  On  the 
abscissa  are  the  possible  values  of  the 
threshold  and  on  the  ordinate  is  the 
fraction  of  points  whose  refinement 
parameter  exceeds  the  selected  threshold. 

There  is  a  knee  in  the  figure  near 
point  C.  For  a  slight  decrease  in  the 
threshold  value,  the  fraction  of  points 
which  would  be  newly  included  in  the 
adaptation  would  rise  rapidly  (a  large 


is  used  to  automatically  set  the  threshold. 


One  must  however  include  guards  into 
this  algorithm.  For  example  it  seems 
foolish  to  choose  a  thresho’*  va  ue  which  is 
less  than  the  average  refi.ie^ent  parameter. 
It  is  also  important  to  ensure  that  too  much 
adaptation  is  not  done  since  that  increases 
the  required  computational  resources.  This 
results  in  constraining  the  possible  values 
of  threshold  by  the  cross-hatched  lines  in 
figure  5. 
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node  number 


Figure  4. 

Distribution  of  refinement  parameter. 


Figure  5. 

Cumulative  distribution  of  refinement 
parameter. 


A  sensitivity  study  was  conducted  to 
ascertain  the  importance  of  the  parameters 
used  for  the  automatic  threshold 
determination.  Figure  6  (a,  b,  and  c)  shows 
the  effect  of  choosing  the  threshold  at 
levels  A,  B,  and  C  respectively.  It  can  be 
seen  that  the  extent  of  the  embedded  region 
Increases  as  the  threshold  value  decreases. 
However,  the  computed  lift  and  drag 
coefficients  show  a  slight  degradation  from 
C  to  B  (about  1  percent)  and  a  larger 
degradation  from  C  to  A  (about  5  percent). 
This  is  consistent  with  the  previous  comment 
that  the  adaptation  is  best  when  it 
completely  surrounds  the  flow  features  of 
Interest. 


a)  Threshold  set  at  A. 


To  determine  the  effectiveness  of  an 
adaptive  grid  scheme,  one  must  consider  its 
accuracy  and  efficiency  as  compared  with 
standard  techniques.  Adaptation 
effectiveness  can  be  measured  in  two  ways; 
the  first  is  with  respect  to  required  CPU 
time  and  the  second  is  with  respect  to 
required  computer  storage.  Adaptation 
accuracy  can  be  measured  by  comparing 
solutions  computed  with  adapted  refinement 
with  those  based  upon  globally  refined 
grids. 

Figure  7  shows  the  computed  accuracy 
versus  efficiency  for  three  different 
airfoils  (to  be  discussed  below).  The 
accuracy  is  measured  by  the  difference 
between  the  computed  lift  coefficient  and  a 
reference  lift  coefficient  (figures  7a  and 
7c)  or  the  difference  between  the  computed 
drag  coefficient  and  a  reference  value 
(figures  7b  and  7d).  The  efficiency  is 
measured  either  by  the  required  CPU  time 
normalized  by  the  time  required  for  the  base 
solution  (figures  7a  and  7b),  or  by  the 
number  of  nodes  similarly  normalized 
(figures  7c  and  7d).  In  each  figure,  for 
each  airfoil  (symbol),  there  are  two  lines. 
The  solid  line  refers  to  adapted  refinement 
and  the  dotted  line  refers  to  global 
refinement. 


The  first  airfoil  is  the  NACA  0012  at 
Mach  number  0.40  and  zero  degrees  angle  of 
attack.  For  this  case,  which  is  denoted  by 
squares,  the  reference  values  of  lift  and 
drag  coefficients  are  both  zero.  It  should 
be  noted  that  all  the  solutions  yielded  zero 
lift,  and  this  is  shown  as  an  error  of  10~ 
(the  assumed  accuracy  of  the  calculation). 
The  drag  coefficients  for  the  solution  with 
one  global  and  one  adapted  refinement  were 
identical,  as  were  the  drags  calculated  for 
two  global  and  two  adapted  refinements. 


b)  Threshold  set  at  B. 


c)  Threshold  set  at  C. 


The  second  case  (denoted  by  circles) 
is  again  for  the  NACA  0012,  but  this  time  at 
Mach  number  0.63  and  2.0  degrees  angle  of 
attack.  For  this  case,  the  reference  lift 
coefficient  is  0.3297  and  the  reference  drag 
coefficient  is  zero.  It  can  again  be  seen 
that  the  lift  and  drag  coefficients  are  the 
same  if  global  or  adapted  refinement  is 


Figure  6.  Effect  of  choice  of  threshold. 

First  difference  of  density. 

RAE  2822,  Mach*0.75,  a*3.0  deg. 


Figure  7.  Comparison  of  accuracy  vs. 
efficiency  for  global  and  adapted  embedding. 


The  third  case  (denoted  by  triangles) 
is  the  RAE  2822  airfoil  at  Mach  number  0.75 
and  3  degrees  angle  of  attack.  The  pressure 
coefficient  distribution  for  this  airfoil 
(figure  8)  shows  a  shock  at  about  75  percent 
chord  on  the  upper  surface,  as  do  the  Mach 
number  contours  for  the  same  case  (figure 
9).  The  three  lines  in  figure  8  are  the 
computed  solutions  without  refinement,  with 
global  refinement,  and  with  adapted 
refinement;  the  global  and  adapted 
refinements  are  so  close  that  they  virtually 
appear  as  one  line  on  the  plot.  For  this 
case,  the  reference  lift  coefficient  is 
1.076  and  the  reference  drag  coefficient  is 
0.0424,  The  anomaly  in  drag  coefficient  for 
the  base  mesh  is  due  to  a  trade-off  between 
shock  location  (and  strength)  change  and 
total  pressure  loss  due  to  smoothing. 
Again,  the  error  in  the  lift  and  drag 
coefficients  are  the  same  for  global  and 
adapted  refinement. 


For  these  cases,  figure  7  demonstrates 
that  adaptive  refinement  yielded  the  same 
accuracy  as  global  refinement,  but  required 
only  12  to  33  percent  as  much  resources. 
Alternatively,  one  notes  that  for  a  given 
resource  allocation,  the  adapted  solutions 
yielded  considerably  more  accurate 
solutions. 


Figure  B.  Pressure  coefficient  for 
non-embedded,  global-embedded,  and 
adapted-embedded  runs. 

RAE  2822  Mach«0.75,  m-3.0  deg. 
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sonic  line 


•  An  adaptive  smoothing  algorithm  has  been 
developed,  which  Is  conservative,  small 
in  smooth  regions,  and  primarily 
diffusive  rather  than  convective. 


•  The  adaptation  strategy  presented  here 
is  extendable  to  three-dimensions,  with 
even  larger  anticipated  savings. 
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The  final  test  case  is  for  the 
supersonic  flow  over  an  8-percent-thick, 
circular-arc,  unstaggered  cascade  at  a  free 
stream  Mach  number  of  1.4.  In  this  case, 
there  Is  a  complex  shock  pattern  which  forms 
from  the  intersection  and  eventual 
coalescence  of  the  leading  and  trailing  edge 
shocks.  It  can  be  seen  that  the  adaptation 
method  used  here  resolves  the  detailed  shock 
interaction  (figure  10a).  If  grid 
refinement  by  redistribution  had  been  used, 
the  case  would  have  yielded  significant 
topological  difficulties.  As  before,  the 
refinement  parameter  is  the  first  difference 
of  density,  with  the  threshold  chosen 
automatically. 


Figure  10b  shows  that  for  the  base 
grid,  the  shock  intersection,  reflection, 
and  eventual  coalescence  is  missed  whereas 
for  either  tne  global  (figure  10c )  or 
adapted  refinement  (figure  lOd),  these 
phenomena  are  obtained  in  detail.  For  this 
case,  the  adapted  solution  ran  4.8  times 
faster  and  required  one  third  the  storage  of 
the  comparable  globally  refined  solution. 

CONCLUSIONS 

•  An  adapted  grid  strategy  which  uses  the 
embedded  mesh  procedure  described  herein 
yields  solutions  with  the  same  accuracy 
as  a  globally  refined  mesh  case  but  only 
requires  between  12  and  33  percent  of 
the  resources  -  for  all  cases  tried  to 
date. 


e  Various  refinement  parameters  have  been 
examined  to  deduce  that  the  first 
difference  of  density  is  the  most 
effective  for  two  dimensional  Euler 
solutions  with  shocks. 


e  A  thresholding  algorithm  has  been 
developed  which  automatically  sets 
thresholds  in  order  to  detect  and 
Isolate  features  from  the  background 
field. 


a)  adapted  grid. 


Base  grid  Time  »  1.00 


c)  Mach  number  contours 
Two  global  refinements.  Time  ■  6.27 


d >  Mach  number  contours 
Two  adapted  refinements.  Time  «  1.31 


Figure  10. 

Solution  for  supersonic  circular  arc  cascade 
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IMPACT  or  SUPERCOMPUTERS  OK  THE  MEET  EEC ACE  CE 
COMPUTATI OKAE  FLUID  DYNAMICS 

Earll  M.  Murman 
Saul  £.  Abarbanel 

Introduction 

A  small  group  of  CFD  researchers  from  the  United  States  and  Israel 
gathered  in  Jerusalem  during  December  1984  at  a  workshop  entitled  “The 
Impact  of  Supercomputers  on  the  Next  Decade  of  Computational  Fluid 
Dynamics."  The  background  of  the  workshop  attendees  ranged  from  CFD 
code  developers  to  applied  mathematicians  to  computer  experts.  During 
the  workshop  the  participants  presented  and  discussed  results  of  their 
current  research.  They  then  encaged  in  discussion  of  the  workshop 
theme.  This  article  attempts  to  summarize  their  observations  and 
speculations  on  what  the  impact  of  supercomputers  will  be  on  CFD  curin' 
the  next  decade.  First,  however,  we  briefly  summarize  the  papers  in 
these  proceedings  and  the  current  status  of  CFD. 

The  Present 

Supercomputers  and  CFD  have  affected  every  aspect  of  fluid 
dynamics  to  some  degree  during  the  past  decade.  Perhaps  the  area 
which  has  experienced  the  most  dramatic  impact  is  the  field  of  attache 
flow  aerodynamics,  typical  of  design  point  conditions  for  transport 
aircraft.  In  this  situation  the  fluid  flow  is  well  behaved  by  design. 
Separated  and  unsteady  flow  are  avoided.  The  turbulent  flow  models 
applicable  to  attached  boundary  layers  are  acceptable  (though  net 
perfect) .  There  are  no  chemical  reactions  or  phase  changes  taking 
place.  The  biggest  challenges  lie  in  solving  the  nonlinear  inviscid 
flow  equations  (primarily  transonic)  and  dealing  with  the  complex 
geometry.  During  the  past  decade,  the  capability  of  the  computers  and 
algorithms  has  developed  enormously  in  this  one  particular  subdisci¬ 
pline.  In  many  instances  they  are  as  much  utilized  as  wind  tunnel 
testing,  although  by  no  means  supplementing  them. 


In  most,  ctner  titles  ct  science  anc  engineering,  many  c:  tne  more 
difficult  fluid  dynamic  phenomena  which  are  absent  in  attained  flow 
aerodynamics  are  of  paramount  importance .  For  example,  turbomachines 
are  doni.nat.ed  X)V  tnre^'Ciinsr.s i onal  viscous  anc  ur.stcacv  phenomena  which 
affect  heat  transfer  and  perforr.ar.ee.  In  many  devices  for  propulsion 
and  chemical  processes,  multicomponent  chemical  reactions  and  turbulent 
mixing  must  be  modeled.  High  performance  aircraft  and  helicopters  are 
strongly  influenced  by  vortical  and  unsteady  effects.  Low  crag  bodies 
are  dominated  by  the  prediction  of  transition.  Separated,  unsteady 
wakes  of  automobiles  influence  their  fuel  consumption  and  handling 
capabilities .  In  large  scale  geophysical  fluid  dynamics,  coriolis 
forces  and  stratification  effects  are  dominant.  These  lead  to  multiple 
time  scale  wave  phenomena.  Unstable  stratification  produces  turbulent, 
buoyant  mixing.  Turbulent  flow  is  present  in  virtually  every'  situation, 
yet  can  only  be  adequately  modeled  for  the  simplest  of  flow’s  like 
attached  boundary  layers  and  jets.  This  does  net  exhaust  the  list,  but 
the  point  made  here  is  quite  clear;  only  the  tip  of  the  iceberg  has 
been  seen  by  the  progress  made  in  attached  flow  aerodynamics.  The 
biggest  challenges  are  yet  to  come.  The  papers  in  the  proceedings  give 
one  assessment  of  where  the  field  stands  in  this  respect. 

The  Pacers 

Fernbach's  paper  gives  a  comprehensive  overview’  of  the  current  per¬ 
formance  of  supercomputers  and  what  is  on  the  horizon.  This  field  is 
now  very  active  following  a  dormant  period  in  the  70’ s.  The  basic 
message  is  that  computer  speed  and  main  memory  will  both  increase  by 
about  two  orders  of  magnitude  in  the  next  decade.  Also,  all  future 
supercomputers  will  be  a  combination  of  vector  (pipeline)  and  parallel 
(multiprocessor)  architectures .  Algorithms  will  have  to  adapt  to  and 
exploit  these  architectural  features  to  achieve  the  stated  machine 
performance.  Thompkins '  paper  demonstrates  that  supercomputing  does 
not  necessarily  have  to  be  done  on  supercomputers .  It  makes  an  inter¬ 
esting  case  for  the  need  cf  personal-sized  supercomputers  with  speeds 
about  one  order  of  magnitude  slower  than  the  mainframes,  but  with  com¬ 
parable  memories.  The  idea  of  using  higher  level  languages  for  multi¬ 
processor  applications  is  introduced  by  Thompkins.  Icavier-Stokes 
results  are  also  presented  for  turbomachine  cascades,  illustrating  the 
state  of  the  art  for  CFD  applied  to  these  flow’s.  The  paper  by  Jameson, 
Leicher,  and  Daw’son  demonstrates  one-  way  that  the  current  generation  cf 


processor  arena recourse. 

The  paper  by  Steger  and  Runzr.g  provides  ar.  overview  cf  current 
issues  regarding  the  computation  cf  irvisric  and  viscous  aerodyr  amic 
flows.  In  addition  to  a  number  cf  interesting  resu_ts  wmcr.  are  pre 
sented,  the  experience  cf  these  authors  using  tne  current  generation 


cf  supercomputers  is  recorded.  Yectoricaticn 


ai  con  trues 


is  explained.  The  need  for  good  graphical  postprocessing  cf  large  data 
bases  turns  out  to  be  mandatory .  P.  similar  message  is  giver,  in  the 
paper  by  Murraan,  Ricci  and  Rowell,  which  compares  tv:o  independently  ob¬ 


tained  solutions  for  leading  edge  vortex  flov 


lor  oelta  wines.  Tr.is 


paper  also  illustrates  that  this  class  cf  compressible  flows  is 
relatively  unexplored  compared  to  shock  wave  dominated  flov:s. 

Several  papers  present  new  algorithms  for  the  Ruler  and  the  incom¬ 
pressible  or  compressible  Kavier-Stokes  equations.  Since  the  solution 
of  these  ecuations  will  become  more  frecuent  with  the  hia'ner  tower 


computers ,  this  is  an  important  topic  for  the  future. 


The  tamer  bv 


Walters  and  Erwoyer  introduces  an  upwind  differencing  line  relaxation 
algorithm  for  the  Euler  equations.  McCormack  presents  a  new  algorithm 
for  the  compressible  Kavier-Stokes  equations  which  has  some  similaritie 
to  the  algorithm  cf  Walters  and  Dwover.  Turkel  presents  methods  for 
accelerating  the  convergence  to  a  steady  state  solution  cf  the  Euler 
or  Kavier-Stokes  equations  by  using  preconditioning  to  alter  the  time 
consistency  of  the  equation  set.  The  paper  by  Glowinski  addresses 
finite  element  methods  for  the  incompressible  Kavier-Stokes  equations 
and  presents  a  number  cf  results  concerning  entry  to  cuots  and  the 
subsequent  internal  flov:.  Remarks  are  also  included  cr.  finite  element 
methods  for  compressible  Kavier-Stokes  equations  and  applications  on 
supercomputers .  Israeli  gives  an  algorithm  for  the  parabolized  Kavier- 
Stokes  (PNS)  equations.  Brandt  and  Ta’asan  present  the  latest  multi¬ 
grid  algorithms  for  quasi-elliptic  systems  which  arise  from  discrete 
approximations  to  the  Kavier-Stokes  and  related  equations.  The 
importance  cf  algorithms ,  such  as  multigrid  methods,  which  have  con¬ 
vergence  rates  (spectral  radius)  independent  cf  the  number  cf  mesh 
points  will  be  mentioned  later. 

Perhaps  no  problem  is  mere  central  to  fluid  mechanics  than  the 
prediction  cf  transition  and  turbulent  flows.  Three  tamers  deal  with 
otic.  Brachet, .  Metcalfe,  Crszag  an 
stability  cf  free  shear  flows  based 
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Che  h a vi e r - S ton e s  equations.  numerical  experiments  such  as  these  car. 
lead  to  new  theoretical  ur.dersoar.d_nq  cf  ir.stahilioy  cf  rotational 
flows.  Feroiger's  paper  gives  a  comprehensive  overview  of  current  and 
future  approaches  to  turbulent  flow  cc reputations  using  direct  simula¬ 
tions  cf  the  havier-Ftok.es  equations,  large  eddy  simulation,  anc 
turbulence  models.  V.’olfshtem  considers  the  latter  topic  ir.  much 
greater  detail.  These  two  papers  point  out  tne  capabilities  and  short¬ 
comings  cf  current  turbulence  models.  The  importance  of  having  accurate 
algorithms  is  stressed  by  both  authors. 

Papers  by  Sulem  and  by  Kichelson  illustrate  how-  numerical  results 
can  be  used  to  understand  the  nature  of  the  solutions  to  partial 
differential  equations.  The  use  of  spectral  methods  for  problems  which 
require  high  accuracy  is  receiving  increased  interest.  Tne  presenta¬ 
tions  by  Abarbanel  and  Gottlieb,  and  Gottlieb  and  Tacmor  consider  some 
basic  issues  regarding  the  resolution  cf  extreme  gradients  by  spectral 
methods.  The  paper  by  Browning  and  Kreiss  illustrates  that  many  fluid 
problems  with  multiple  time  and  length  scales  are  exceedingly  difficult 
to  compute,  even  with  "unlimited"  computer  power.  It  is  important  to 
understand  that  the  powerful  new  supercomputers  will  only  yield  useful 
results  if  the  mathematical  and  numerical  analysis  formulation  is  care¬ 
fully  done.  The  paper  by  Sever  is  another  illustration  of  this. 


The  Next  Decade 

Turing  the  next  decade  supercomputer  power  will  increase  dramatical¬ 
ly.  The  directly  addressable  high  speed  memory  capacity  will  increase 
by  about  two  orders  cf  magnitude,  from  2  million  words  to  256  million 
words.  Processor  speed  will  increase  an  order  cf  magnitude  from  about 
100  MFLOPS  to  1000  KFL0PS,  or  perhaps  more.  It  is  likely  that  the 
corresponding  parameters  cf  smaller  computers  will  increase  by  similar 
factors.  History  has  shown  that  whenever  an  important  parameter  is 
varied  by  an  order  of  magnitude,  new  discoveries  are  made.  The  diffi¬ 
culty  is  to  have  some  feeling  as  to  what  these  discoveries  might  be. 

Tne  participants  realized,  cf  course,  that  forecasting  the  future  is 
never  an  exact  science,  and  one  should  proceed  cautiously.  it  is 
interesting  to  note,  however,  that  during  t he  panel  discussion  almost 
everyone  subscribed  to  the  idea  that  the  new  supercomputers  will  net 
or.Iv  allow  tackling  bigger  problems,  but  will  also  lead  to  a  better 
understanding  cf  the  physics  of  some  complex  problems  such  as  turbulence, 
vertical  flows ,  and  chemically  reacting  flows . 


.vv'y.,y..v.yn  {•  ^ 


_e  v:e  Eu.~a 


enaees  co 


'jest  ions  v:r. 


Suoerccn.r 


In  tne  tiexc  c:  aerocynanrj.es,  one  p reiir.cncry  aestgr.  c:  transport 
aircraft  will  primarily  be  cone  or.  supercomputers.  The  modeling  anc 
computing  capability  will  be  basically  in  place.  Unlike  earlier  pre¬ 
dictions  that  computers  will  make  wind  tunnels  obsolete,  few  people 
subscribe  to  that  viewpoint  now.  What  is  more  likely  to  happen  is  that 
the  use  of  wind  tunnels  by  researchers  and  aesicn  encineers  will  cnanae 


Less  and  less  of  the  exploratory  design  will  be  done  by  tests  as  the 
oredictive  methods  become  more  reliable.  This  has  already  natoened  in 


several  instances  with  the  current  aeneration  of  commuters. 


The  next 


generation  will  provide  enough  resolution  and  speed  that  a  realistic 
model  of  an  actual  cruising  transport  aircraft  can  be  computed. 

The  capability  to  model  "off-design"  or  "unclean"  aerodynamic  flows 
will  increase.  These  are  flows  which  are  separated,  unsteady,  vortex 
dominated,  and  the  like.  Such  flow’s  are  cf  great  importance  for 
maneuvers  of  hich  oerfcrtiance  aircraft  or  for  emeraency  situations  for 


iranscort  aircraft 


The  loads  aevelooed  in  these  recimes  often 


determine  tne  required  strength  cr  the  aircraft  components.  The  same 
phenomena  often  dominate  rotary’  wing  and  rotating  machinery’  aero¬ 
dynamics  .  Caoacitv  cf  computers  and  alccrithms  ur>  to  now  has  not  been 


aaecuate  to  suet 


frontal  assault  on  this  class  cf  crchlems.  The 


payoff  for  analy’sis  cf  unsteady’,  separated,  vortical  flows  will  be  much 
greater  than  for  the  clean  flows  representing  cruise  aerodynamics. 

This  is  because  little  theory’  has  ever  been  developed  for  them.. 

The  complexity’  cf  problems  which  the  researcher  and  the  engineer 
v.’ill  be  dealing  with  will  crow’  in  some  proportion  to  the  new  computer 
power.  This  will  have  a  number  of  impacts  cn  the  daily  life  cf  the 
fluids  mechanician.  Problems  'under  investigation  will  have  many’  length 
and  time  scales.  Analy’sis  cf  results  will  be  more  challenging.  It 
will  be  harder  to  understand  the  solution  due  to  the  number  cf  inter¬ 
acting  phy’sical  phenomena  present.  This  is  illustrated  by’  the  paper 
cf  Eraohet,  Metcalfe,  Orstag,  and  Riley.  Kew  methods  cf  analyzing  and 
presenting  results  will  have  to  eneroe  in  order  to  deal  with  this. 


meat  output  is  cruet 
ethnology  will  nelo  o 


Cue  difficulty  which  can  be  fcreseer.  is  the  problem  cf  verifying 
the  accuracy  or  fidelity  cf  a  computed  result.  Up  to  now,  it  has 
generally  been  possible  to  compare  computed  results  with  theory  for 
limiting  conditions.  For  example,  a  transonic  v:ing  calculation  can  be 
compared  with  linear  wing  theory  for  low  Mach  numbers,  or  a  Havier- 
Stokes  solution  can  be  compared  with  a  laminar  boundary  layer.  But  as 
the  computations  move  into  more  nonlinear  flows,  the  past  theoretical 
framework  will  become  less  and  less  applicable.  Comparison  v.*ith  ex¬ 
periment  is  an  essential,  and  independent  computations  cf  the  same 
problem  by  different  researchers  will  be  necessary.  Perhaps  a  renewed 
interest  in  theory  will  result  from  this  need. 

Imoact  of  Supercomputers  on  Basic  Sciences 

As  one  participant  stated,  the  great  masters  of  fluid  mechanics  in 
the  past  solved  all  the  linear  problems  and  left  us  with  only  the  non¬ 
linear  ones.  Since  most  fluid  mechanic  problems  are  nonlinear,  we  can 
speculate  that  the  ability  to  model  highly  nonlinear  problems  with 
powerful  computers  will  lead  to  many  new  discoveries.  Another  partici¬ 
pant  thought  that  the  impact  cf  supercomputers  will  influence  the  basic 
way  we  think  about  physical  problems.  Hew  information  will  be  dis¬ 
covered  from  numerical  experiments  and  provide  insight  for  modeling. 

In  this  sense,  computational  experiments  are  akin  to  laboratory  experi¬ 
ments  which  have  provided  insight  and  ideas  throughout  the  history’  cf 
fluid  mechanics  and  ether  scientific  disciplines. 

In  the  past,  computational  methods  have  made  a  major  impact  on  our 
ability  to  compute  and  understand  potential  flows  and  inviscia  flow’s 
dominated  by  shockwaves.  One  can  conclude  that  the  classes  of  flows 
are  well  understood  both  from  the  physical  and  algorithmic  points  of 
view.  Although  the  ability  to  analyze  shock  dominated  flow’s  has  been  a 
major  step  forward  in  fluid  mechanics,  much  is  left  to  be  acne.  For 
example,  only  limited  CTD  studies  have  been  done  for  vertical  flows, 
and  I'tile  is  understood  about  the  algorithm  requirements  for  inviscid 
rotational  flow’s.  Many  studies  have  been  done  for  two-dimensional 
separated  flow’s,  but  only  limited  studies  for  three-dimensional 
separated  flow’s.  Although  efficient  algorithms  for  steady  flows  are 
under  development,  indications  are  that  the  flows  these  algorithms  are 
to  be  applied  to  may  be  unsteady  in  nature.  See  fur  example  Thompkins 
or  Murman,  Rizzi  and  Powell. 

Perhaps  no  area  is  more  tempting  to  speculate  on  than  the  field  cf 
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nant  since  Reynolds  introduced  all  the  unknowns  without  introducing  any 
new  equations.  In  the  past  decade,  computational  models  and  laboratory 
experiments  have  opened  a  new  look  at  turbulence.  The  idea  of  organized 
or  coherent  structure  has  emerged.  On  the  ether  hand,  mathematicians 
have  shown  that  solutions  to  fairly  simple  dynamical  systems  have 
chaotic  behavior.  An  interesting  question  which  v;as  posed  is  "A’hat 
will  be  the  resolution  of  the  speculation  that  there  is  both  determinism 
and  chaos  in  nonlinear  equations?"  Computational  experiments  could 
provide  a  framework  for  helping  to  answer  this  question. 

Another  area  which  will  probably  be  strongly  influenced  by  mere 
powerful  computational  approaches  is  the  coupling  of  chemical  reactions 
and  heat  release  to  fluid  flow  problems.  Even  fairly  "simple"  reactions 
involve  many  species  with  many  rime  scales  of  reactions.  In  the  past, 
computers  simply  were  not  large  enough  to  tackle  many  of  these  problems. 
Rate  constants  are  always  an  uncertain  factor  in  such  calculations. 
Perhaps  being  able  to  model  the  experimental  conditions  under  w’hich  the 
rate  constants  are  measured  will  lead  to  more  accurate  measurements  of 
their  values. 

One  issue  on  which  there  was  quite  a  difference  of  opinion  is  the 
degree  to  which  modeling  will  be  required  prior  to  computing.  On  the 
one  hand,  many  participants  felt  that  the  time  was  upon  us  to  tackle 
the  full  three-dimensional  IC  avis  r- Stokes  equations,  possibly  adding 
models  only  fer  subgrid  scale  turbulence.  Others  felt  that  the  past 
practice  of  selecting  simplified  sets  cf  equations  such  as  inviscia  or 
parabolized  viscous  will  still  be  prudent.  It  is  likely  that  some 
level  of  modeling  will  always  be  required  as  computer  capability  will 
never  be  big  enough  to  solve  a  complex  problem  from  first  principles. 

In  fact,  for  most  problems  this  is  unnecessary.  The  question  is, 
will  the  type  of  modeling  appropriate  for  the  future  be  different  from 
that  used  in  the  past  when  computer  memory,  speed,  and  accessibility 
were  more  limited? 

Imsact  of  Supercomputers  on  Alccrithms  and  Languages 

An  important  issue  regarding  algorithms  arises  from  the  multi¬ 
processor  and  vector  architectures  cf  supercomputers .  Algorithms  which 
cannot  be  efficiently  used  or.  these  architectures  will  be  cf  limited 
utility,  ti any  fluid  mechanic  problems  are  solved  using  time  dependent 
integration  procedures  for  initial  boundarv  value  problems.  Both 
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explicit,  and  implicit  methods  are  utilized.  lr.  ceneral,  explici-t 
algorithms  are  easier  to  vecocrrze  than  tr.p licit  ones.  The  latter 
usually  involve  recursive  steps  ii.  the  matrix  inversion.  However,  the 
paper  bv  Steger  and  Buning  demonstrates  an  effective  vectcrization 
strategy  for  simultaneous  inversion  cf  a  large  number  cf  tridiagonal 
matrices  applicable  to  approximate  factorization  methods.  Explicit 
algorithms  are  also  easier  to  adapt  to  multiprocessor  architectures  as 
the  solution  domain  can  be  subdivided  without  introducing  difficulties 
in  the  algorithm.  The  paper  by  Jameson,  Leicher  and  Dawson  reports  a 
strategy  and  results  which  illustrate  this.  Effective  strategies  for 
adapting  implicit  algorithms  for  multiprocessor  architectures  need  to 
be  found  also.  Participants  generally  agreed  that  the  real  payoff  is 
for  algorithms  which  can  work  effectively  on  tens  or  hundreds  of 
parallel  processors,  not  just  two,  four,  or  eight. 

It  is  clear  that  algorithm  developers  must  be  cognizant  of  the 
advantages  and  constraints  which  non-Von  Neumann  architectures  will 
place  on  supercomputing.  The  paper  by  Thompkins  is  an  indicator  of 
what  will  come.  In  addition  to  the  conceptual  changes  in  algorithms 
due  to  multiprocessors  and  vector  processors,  efficiency  limitations 
arise  in  the  speed  at  which  main  memory  can  be  accessed  from  various 
processors  or  the  speed  at  which  data  can  be  exchanged  between  pro¬ 
cessors.  The  personal-size  supercomputer  reported  by  Thompkins  involves 
two  processors  (a  scalar  host  and  an  attached  vector).  Movement  cf  the 
data  must  be  carefully  managed  to  avoid  bottlenecks.  In  the  future 
algorithm  development  must  take  into  consideration  not  only  traditional 
numerical  analysis,  but  also  computer  science  aspects. 

Another  algorithm  issue  arises  from  the  shear  size  cf  the  main 
memory'  cf  supercomputers .  Problems  with  very  fine  meshes  will  be 
possible.  Kith  the  exception  cf  multigrid  algorithms  for  elliptic 
equations,  the  asy'mptotic  convergence  rate  (spectral  radius)  of 
iterative  methods  is  dependent  on  the  number  cf  mesh  points.  Recent 
estimations  done  by'  researchers  at  KA.SA  Eangiey'  indicate  that  the  time 
required  to  reach  convergence  for  three-dimensional  Navier-Srokes 
equations  on  a  256  megaword  machine  is  excessive  to  the  point  cf  being 
unrealistic.  This  indicates  a  real  need  for  finding  algorithms  which 
have  asymptotic  convergence  rates  which  are  either  mesh  independent  or 
vary  (at  worst)  slowly  with  the  number  cf  grid  points.  Such  algorithms 
must  work  for  problems  with  widely'  varying  length  and  time  scales,  not 
just  model  problems  cr  model  equations  (see  Jameson  et  al) . 
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Tnere  was  considerable  discussror.  or.  tne  need  .for  hi  oner  level 
languages  mat  will  make  it  easier  to  construct  a  scluttor.  accroach  for 
a  new  problem,  as  well  as  make  it  easier  to  utilise  the  new  architecture 
The  general  strategy  for  solving  a  pretier,  ry  IFF  as  more  cr  less 
common,  A  grid  is  generated,  discretisation  cf  spatial  derivatives  and 
boundary  conditions  is  done,  an  iteration  cr  time  integration  method  is 
selected,  and  various  outputs  are  required.  It  was  suggested  that 
assembling  these  tools  and  manipulating  them  for  embedded  subdomains 
and  the  like  would  lend  itself  to  a  higher,  and  therefore  simpler, 
language.  FORTRAN  is  the  language  cf  the  CFD  community  to  date.  The 
paper  by  Tnompkins  introduces  some  higher  level  constructs  for  managing 
the  solution  process  in  a  multiprocessor  environment. 

Another  area  in  which  algorithm  innovation  may  be  required  is  in 
preprocessing  (grid  generation)  and  post-processinq  (data  base  analysis) 
The  papers  by  Steger  and  Buning  anc  Murman,  Ricci  and  Pow’ell  indicate 
that'  graphical  analysis  is  imperative,  but  other  ways  cf  manipulating 
the  data  bases  would  be  desirable.  Maybe  knowledge  base  programs 
("expert  systems")  will  be  helpful  in  finding  the  important  features  in 
a  given  solution.  Or  perhaps  pattern  recognition  approaches  will  be 
needed. 

Impact  on  Subsystems 

The  workshop  attendees  for  the  most  part  represented  users  cf 
supercomputers ,  and  net  hardware  specialists.  However,  with  the  large 
data  bases  which  will  be  generated,  the  participants  felt  that  several 
cf  the  supporting  subsystems  might  vjell  be  inadequate  to  match  tne 
power  cf  the  high  speed  processors.  Participants  who  have  had  experi¬ 
ence  with  supercomputers  felt  that  strong  graphics  capability  is  a 
number  one  requirement  for  analycing  results.  As  discussed  above, 
some  new  graphics  algorithms  may  well  need  to  be  devised  to  analyse 
complicated  flow  fields.  But  powerful  processing  and  graphical  display 
capabilities  are  also  required.  As  the  paper  by  Tnompkins  points  out, 
a  researcher  will  typically  spend  much  cf  his  cr  her  time  performing 
graphical  analysis  which  requires  a  processor  about  an  order  cf  mag¬ 
nitude  smaller  than  the  supercomputer .  A  large  number  cf  operations 
are  required  to  construct  the  output  quantities  which  typically  are 
different  combinations  of  the  dependent  variables  than  those  which  are 
stored  in  the  data  base.  Many’  researchers  currently’  thin):  that  t 
plot  files  will  be  created  on  the  supercomputer  itself  for  these 
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reasons.  However,  super  graphical  processor.-, ,  which  coulc  be  much 
cheaper  and  therefore  rr.ore  available  ic  the  users,  should  be  develops: 
Typically ,  algorithms  for  creating  the  graphical  data  base  are  easily 
vectorized  and  could  be  done  on  array  processors  attached  to  most  pro 
cessors.  Tnese  could  lead  to  very  powerful  and  inexpensive  graphics 
workstations.  Fortunately  the  technology  is  developing  rapidly  for 
medium  to  high  resolution  color  display  devices  v:ith  interactive 

capability. 

Most  of  the  attendees  at  the  workshop  fall  in  the  category  of 
"remote"  users.  They  are  not  located  at  the  same  site  as  the  super¬ 
computer.  There  was  significant  concern  that  remote  communications 
will  be  a  real  bottleneck.  Regular  dial-up  capability  is  barely 
adequate  at  present  for  editing  files  or  transmitting  small  output 
files  to  remote  users.  It  certainly  would  be  impossible  to  do  inter¬ 
active  graphics  processing  from  a  remote  site  or  to  transmit  the 
entire  data  base  -for  onsite  analysis  using  even  dedicated  data  lines 
currently  available.  The  best  w_y  to  communicate  with  remote  facili¬ 
ties  at  present  is  to  transmit  magnetic  tapes  via  express  mail 
services.  This  inevitably  leads  to  delays  and  slow  turnaround.  The 
impact  of  supercomputers  on  researchers  who  are  not  co-locatea  with 
the  machines  will  be  minor  if  high  bandwidth  communications  are  not 
available. 

Conclusions 

The  impact  of  supercomputers  on  the  next  decade  cf  computational 
fluid  dynamics  will  be  substantial.  Kith  processor  speeds  and  high 
speed  memory  increasing  by  two  orders  of  magnitude,  many  changes  will 
take  place.  The  difficulty  lies  in  accurately  forecasting  what  those 
changes  will  be.  The  above  discussion  presents  the  thinking  cf  one 
group  of  active  CTD  researchers.  Perhaps  their  viewpoints  will  sene 
to  help  others  to  become  aware  of,  and  think  about,  these  changes  as 
they  take  place.  The  next  decade  has  already  started, 
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